include("Flux.jl")

function EulerRhs(Q,N,K,vmapM,vmapP,mapI, mapO, mapW,mapC,Fx,Fy,Drw,Dsw,rx,sx,ry,sy,nx,ny,time,gamma,LIFT,Fscale)
    """
    Purpose: Evaluate RHS in 2D Euler equations, discretized on weak form
    with a numerical flux
    """
    Nfaces=3; Nfp = N+1; Np = convert(Int, (N+1)*(N+2)/2)
    vmapM = reshape(vmapM, Nfp*Nfaces, K)
    vmapP = reshape(vmapP, Nfp*Nfaces, K)

    # Compute volume contributions
    F,G,rho,u,v,p = EulerFluxes(Q, gamma)

    # Compute weak derivatives
    rhsQ = zeros(Float64, Np,K,4)
    for n=1:4
        dFdr = Drw*F[:,:,n]; dFds = Dsw*F[:,:,n]
        dGdr = Drw*G[:,:,n]; dGds = Dsw*G[:,:,n]
        rhsQ[:,:,n] = (rx.*dFdr + sx.*dFds) + (ry.*dGdr + sy.*dGds)
    end

    # Evaluate '-' and '+' traces of conservative variables
    QM = zeros(Float64,size(vmapM)[1],size(vmapM)[2],4)
    QP = zeros(Float64,size(vmapP)[1],size(vmapP)[2],4)
    for n=1:4
        Qn = Q[:,:,n]
        QM[:,:,n] = Qn[vmapM]; QP[:,:,n] = Qn[vmapP]
    end
    # set boundary conditions by modifying positive traces
    QP = ExactSolutionBC(Q, Np, K, Fx, Fy, nx, ny, mapI, mapO, mapW, mapC, QP, time)

    # Evaluate primitive variables & flux functions at '-' and '+' traces
    fM,gM,rhoM,uM,vM,pM = EulerFluxes(QM, gamma)
    fP,gP,rhoP,uP,vP,pP = EulerFluxes(QP, gamma)

    # Compute local Lax-Friedrichs/Rusonov numerical fluxes
    lambda = max.(sqrt.(uM.^2+vM.^2) + sqrt.(abs.(gamma*pM./rhoM)),
                  sqrt.(uP.^2+vP.^2) + sqrt.(abs.(gamma*pP./rhoP)))
    lambda = reshape(lambda, Nfp, Nfaces*K)
    lambda = ones(Nfp, 1)*maximum(lambda, dims = 1)
    lambda = reshape(lambda, Nfp*Nfaces, K)

    # Lifting fluxes
    for n=1:4
      nflux = nx.*(fP[:,:,n] + fM[:,:,n]) + ny.*(gP[:,:,n] + gM[:,:,n]) +
              lambda.*(QM[:,:,n] - QP[:,:,n])
      rhsQ[:,:,n] = rhsQ[:,:,n] - LIFT*(Fscale.*nflux/2)
    end
    rhsQ
end
